# Loughney et al., "Tectonic influence on Cenozoic mammal richness and sedimentation history of the Basin and Range, western North America"
# Science Advances
# R code to plot Figure 4 A, B, and C
# final published figures processed in Adobe Illustrator CC

# load data files
macrostrat <- read.csv("macrostrat_June-2021.csv", header = TRUE, stringsAsFactors = FALSE)
fossils <- read.csv("fossiliferous_units_January-2021.csv", header = TRUE, stringsAsFactors = FALSE)
area <- read.csv("area_change.csv", header = TRUE, stringsAsFactors = FALSE)
deformation <- read.csv("deformation_rates_January-2021.csv", header = TRUE, stringsAsFactors = FALSE)

baseAge <- 36
topAge <- 0

# ------------------------------------------------------------------------------------------------------------------------
# Figure 4 B - deformation rates of subregions of the Basin and Range
# plot
plot.new();
par(mfrow = c(1,1))
par(oma = c(1, 1, 1, 1));
par(mar = c(5, 4, 4, 2));
plot.window(xlim = c(37, -0.5), ylim = c(-30, 190));
axis(side = 1, at = seq(0, 36, by = 2), cex.axis = 0.9, pos = -16)
mtext("Age (Ma)", side = 1, line = 1, cex.lab = 1)
axis(side = 2, at = seq(0, 190, by = 10), cex.axis = 0.9, pos = 37, las = 1)
mtext("Total deformation rate (km/Myr)", side = 2, line = 3, cex.axis = 1)
# time scale
rect(xleft = c(baseAge, baseAge), ybottom = c(-15, -15), xright = c(33.9, 33.9), ytop = c(-5, -5)) # Eocene
text(x = baseAge-(baseAge-33.9)/2, y = -10, "Eo", cex = 0.9)
rect(xleft = c(33.9, 33.9), ybottom = c(-15, -15), xright = c(23.03, 23.03), ytop = c(-5, -5)) # Oligocene
text(x = 33.9-(33.9-23.03)/2, y = -10, "Oligocene", cex = 0.9)
rect(xleft = c(23.03, 23.03), ybottom = c(-15, -15), xright = c(5.33, 5.33), ytop = c(-5, -5)) # Miocene
text(x = 23.03-(23.03-5.33)/2, y = -10, "Miocene", cex = 0.9)
rect(xleft = c(5.33, 5.33), ybottom = c(-15, -15), xright = c(1.8, 1.8), ytop = c(-5, -5)) # Pliocene
text(x = 5.33-(5.33-1.8)/2, y = -10, "Plio", cex = 0.9)
rect(xleft = c(1.8, 1.8), ybottom = c(-15, -15), xright = c(0, 0), ytop = c(-5, -5)) # Pleistocene
text(x = 1.8-(1.8-topAge)/2, y = -10, "Q", cex = 0.9)
# deformation curves
lines(deformation$mid_bin, deformation$NB_rate, col = "#65A9DD", lwd = 2, lend = "butt")
polygon(x = c(deformation$mid_bin, rev(deformation$mid_bin)), y = c(deformation$NB_upper, rev(deformation$NB_lower)), col = adjustcolor("#65A9DD", alpha.f = 0.30), border = NA)
lines(deformation$mid_bin, deformation$CB_rate, col = "#CD4F39", lwd = 2, lend = "butt")
polygon(x = c(deformation$mid_bin, rev(deformation$mid_bin)), y = c(deformation$CB_upper, rev(deformation$CB_lower)), col = adjustcolor("#CD4F39", alpha.f = 0.30), border = NA)
lines(deformation$mid_bin, deformation$SB_rate, col = "#FEC124", lwd = 2, lend = "butt")
polygon(x = c(deformation$mid_bin, rev(deformation$mid_bin)), y = c(deformation$SB_upper, rev(deformation$SB_lower)), col = adjustcolor("#FEC124", alpha.f = 0.30), border = NA)
legend(34, 190, c("Northern BR", "Central BR", "Southern BR"), c("#65A9DD", "#CD4F39", "#FEC124"))

# ----------------------------------------------------------------------------------------------------------------------
# Figure 4 C - sediment-accumulation rates of fossiliferous units of subregions of the Basin and Range
plot.new();
par(mfrow = c(1,1))
par(oma = c(1, 1, 1, 1));
par(mar = c(5, 4, 4, 2));
plot.window(xlim = c(37, -0.5), ylim = c(-500, 3000));
axis(side = 1, at = seq(0, 36, by = 2), cex.axis = 0.9, pos = -260)
mtext("Age (Ma)", side = 1, line = 1, cex.lab = 1)
axis(side = 2, at = seq(0, 3000, by = 500), cex.axis = 0.9, pos = 37, las = 1)
mtext("Sediment-accumulation rate (m/Myr)", side = 2, line = 3, cex.axis = 1)
# time scale
rect(xleft = c(baseAge, baseAge), ybottom = c(-250, -250), xright = c(33.9, 33.9), ytop = c(-75, -75)) # Eocene
text(x = baseAge-(baseAge-33.9)/2, y = -150, "Eo", cex = 0.9)
rect(xleft = c(33.9, 33.9), ybottom = c(-250, -250), xright = c(23.03, 23.03), ytop = c(-75, -75)) # Oligocene
text(x = 33.9-(33.9-23.03)/2, y = -150, "Oligocene", cex = 0.9)
rect(xleft = c(23.03, 23.03), ybottom = c(-250, -250), xright = c(5.33, 5.33), ytop = c(-75, -75)) # Miocene
text(x = 23.03-(23.03-5.33)/2, y = -150, "Miocene", cex = 0.9)
rect(xleft = c(5.33, 5.33), ybottom = c(-250, -250), xright = c(1.8, 1.8), ytop = c(-75, -75)) # Pliocene
text(x = 5.33-(5.33-1.8)/2, y = -150, "Plio", cex = 0.9)
rect(xleft = c(1.8, 1.8), ybottom = c(-250, -250), xright = c(0, 0), ytop = c(-75, -75)) # Pleistocene
text(x = 1.8-(1.8-topAge)/2, y = -150, "Q", cex = 0.9)
# SAR
lines(fossils$mid_bin, fossils$NB_SAR, type = "S", col = "#65A9DD", lwd = 2, lend = "butt")
lines(fossils$mid_bin, fossils$CB_SAR, type = "S", col = "#CD4F39", lwd = 2, lend = "butt")
lines(fossils$mid_bin, fossils$SB_SAR, type = "S", col = "#FEC124", lwd = 2, lend = "butt")
legend(34, 3000, c("Northern BR", "Central BR", "Southern BR"), c("#65A9DD", "#CD4F39", "#FEC124"))

# ----------------------------------------------------------------------------------------------------------------------
# Figure 4 D - area-change rates of subregions of the Basin and Range
plot.new();
par(mfrow = c(1,1))
par(oma = c(1, 1, 1, 1));
par(mar = c(5, 4, 4, 2));
plot.window(xlim =c(37, -0.5), ylim = c(-400, 2000));
axis(side = 1, at = seq(0, 36, by = 2), cex.axis = 0.9, pos = -210)
mtext("Age (Ma)", side = 1, line = 1, cex.lab = 1)
axis(side = 2, at = seq(0, 2000, by = 500), cex.axis = 0.9, pos = 37, las = 1)
mtext(expression("Area-change rate (km"^2*"/Myr)"), side = 2, line = 3, cex.axis = 1)
# time scale
rect(xleft = c(baseAge, baseAge), ybottom = c(-200, -200), xright = c(33.9, 33.9), ytop = c(-60, -60)) # Eocene
text(x = baseAge-(baseAge-33.9)/2, y = -130, "Eo", cex = 0.9)
rect(xleft = c(33.9, 33.9), ybottom = c(-200, -200), xright = c(23.03, 23.03), ytop = c(-60, -60)) # Oligocene
text(x = 33.9-(33.9-23.03)/2, y = -130, "Oligocene", cex = 0.9)
rect(xleft = c(23.03, 23.03), ybottom = c(-200, -200), xright = c(5.33, 5.33), ytop = c(-60, -60)) # Miocene
text(x = 23.03-(23.03-5.33)/2, y = -130, "Miocene", cex = 0.9)
rect(xleft = c(5.33, 5.33), ybottom = c(-200, -200), xright = c(1.8, 1.8), ytop = c(-60, -60)) # Pliocene
text(x = 5.33-(5.33-1.8)/2, y = -130, "Plio", cex = 0.9)
rect(xleft = c(1.8, 1.8), ybottom = c(-200, -200), xright = c(0, 0), ytop = c(-60, -60)) # Pleistocene
text(x = 1.8-(1.8-topAge)/2, y = -130, "Q", cex = 0.9)
# area curves
lines(deformation$mid_bin, area$NB_rate_change_sum, col = "#65A9DD", lwd = 2, lend = "butt")
lines(deformation$mid_bin, area$CB_rate_change_sum, col = "#CD4F39", lwd = 2, lend = "butt")
lines(deformation$mid_bin, area$SB_rate_change_sum, col = "#FEC124", lwd = 2, lend = "butt")
legend(34, 2000, c("Northern BR", "Central BR", "Southern BR"), c("#65A9DD", "#CD4F39", "#FEC124"))

